#ifndef _MATRIX_H_
#define _MATRIX_H_

#include <iostream>

template <class T>
class TMatrix {
private:
  int nn;                       // column
  int mm;                       // row
  T **v;
public:
  TMatrix();
  TMatrix(int m, int n);
  TMatrix(int m, int n, const T &a);
  TMatrix(int m, int n, const T *a);
  TMatrix(const TMatrix &rhs);

  TMatrix& operator=(const TMatrix &rhs);
  inline T* operator[](const int i);
  inline const T* operator[](const int i) const;
  inline int nrows() const;
  inline int ncols() const;
  void resize(int newm, int newn);
  void assign(int newm, int newn, const T &a);
  void print(void); 
  ~TMatrix();
};

template <class T>
TMatrix<T>::TMatrix():mm(0), nn(0), v(NULL) {}

template <class T>
TMatrix<T>::TMatrix(int m, int n)
  : mm(m), nn(n), v(m > 0?new T*[m]:NULL)
{
  int i, nel = m*n;
  if (v) v[0] = nel > 0? new T[nel]:NULL;
  for (i = 1; i < m; i++)
    v[i] = v[i-1] + n;
}

template <class T>
TMatrix<T>::TMatrix(int m, int n, const T &a)
  :mm(m), nn(n), v(m>0?new T*[m]:NULL)
{
  int i, j, nel = m*n;
  if (v) v[0] = nel > 0 ? new T[nel] : NULL;
  for (i = 1; i < m; i++)
    v[i] = v[i-1] + n;
  for (j = 0; j < m; ++j)
    for (i = 0; i < n; ++i)
      v[j][i] = a;
}

template <class T>
TMatrix<T>::TMatrix(int m, int n, const T *a)
  :mm(m), nn(n), v(m>0?new T*[m]:NULL)
{
  int i, j, nel = m*n;
  if (v) v[0] = nel > 0 ? new T[nel]:NULL;
  for (i = 1; i < m; i++)
    v[i] = v[i-1] + n;
  for (j = 0; j < m; j++)
    for (i = 0; i < n; i++)
      v[j][i] = *a++;
}

template <class T>
TMatrix<T>::TMatrix(const TMatrix &rhs)
  :mm(rhs.mm), nn(rhs.nn), v(mm>0?new T*[mm]:NULL)
{
  int i, j, nel = mm*nn;
  if (v) v[0] = nel>0?new T[nel]:NULL;
  for (i=1; i < mm; i++)
    v[i] = v[i-1] + nn;
  for (j = 0; j < mm; j++)
    for (i = 0; i < nn; i++)
      v[j][i] = rhs[j][i];
}

template <class T>
TMatrix<T>& TMatrix<T>::operator=(const TMatrix<T> &rhs)
{
  if (this != &rhs) {
    int i, j, nel;
    if (nn != rhs.nn || mm != rhs.mm) {
      if (v != NULL) {
        delete [] (v[0]);
        delete [] (v);
      }
      nn = rhs.nn;
      mm = rhs.mm;
      v = mm > 0 ? new T*[mm] : NULL;
      nel = mm * nn;
      if (v) v[0] = nel > 0 ? new T[nel] : NULL;
      for (i = 1; i < mm; i++)
        v[i] = v[i-1] + nn;
    }
    for (j = 0; j < mm; j++)
      for (i = 0; i < nn; i++)
        v[j][i] = rhs[j][i];
  }
  return *this;
}

template <class T>
inline T* TMatrix<T>::operator[](const int i)
{
  return v[i];
}

template <class T>
inline const T* TMatrix<T>::operator[](const int i) const
{
  return v[i];
}

template <class T>
inline int TMatrix<T>::nrows() const
{
  return mm;
}

template <class T>
inline int TMatrix<T>::ncols() const
{
  return nn;
}

template <class T>
void TMatrix<T>::resize(int newm, int newn)
{
  int i, nel;
  if (newm != mm || newn != nn) {
    if (v != NULL) {
      delete [] (v[0]);
      delete [] (v);
    }
    mm = newm;
    nn = newn;
    v = mm > 0? new T*[mm] : NULL;
    nel = mm*nn;
    if (v) v[0] = nel>0?new T[nel]:NULL;
    for (i=1; i <mm; i++)
      v[i] = v[i-1] + nn;
  }
}

template <class T>
void TMatrix<T>::assign(int newm, int newn, const T &a)
{
  int i, j, nel;
  resize(newm, newn);
  for (j = 0; j < mm; j++)
    for (i = 0; i < nn; i++)
      v[j][i] = a;
}

template <class T>
TMatrix<T>::~TMatrix()
{
  if (v != NULL) {
    delete [] (v[0]);
    delete [] (v);
  }
}

template <class T>
void TMatrix<T>::print(void)
{
  int i, j;
  for (j = 0; j < mm; j++) {
    for (i = 0; i < nn; i++)
      std::cout << v[j][i] << " ";
    std::cout << std::endl;
  }
}

typedef TMatrix<double> MatDoub;
typedef TMatrix<int> MatInt;

#endif /* _MATRIX_H_ */
